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Abstract 

We study the MSSM finite temperature electroweak phase transition with lattice Monte 
Carlo simulations, for a large Higgs mass (uih ~ 95 GeV) and light stop masses 
( m tR ~ . . . 160 GeV). We employ a 3d effective field theory approach, where the 
degrees of freedom appearing in the action are the SU(2) and SU(3) gauge fields, 
the weakly interacting Higgs doublet, and the strongly interacting stop triplet. We 
determine the phase diagram, the critical temperatures, the scalar field expectation 
values, the latent heat, the interface tension and the correlation lengths at the phase 
transition points. Extrapolating the results to the infinite volume and continuum limits, 
we find that the transition is stronger than indicated by 2-loop perturbation theory, 
guaranteeing that the MSSM phase transition is strong enough for baryogenesis in this 
regime. We also study the possibility of a two-stage phase transition, in which the 
stop field gets an expectation value in an intermediate phase. We find that a two-stage 
transition exists non-perturbatively, as well, but for somewhat smaller stop masses than 
in perturbation theory. Finally, the latter stage of the two-stage transition is found to 
be extremely strong, and thus it might not be allowed in the cosmo logical environment. 
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1 Introduction 



The electroweak phase transition is the last instance in the history of the Universe 
that a baryon asymmetry could have been generated |], and as such, also the scenario 
requiring the least assumptions beyond established physics. In principle, even the 
Standard Model contains the necessary ingredients for baryon number generation (for 
a review, see ||). However, on a more quantitative level, the Standard Model is too 
restricted for this purpose: for the allowed Higgs masses m H > 75 GeV, it turns out 
that there would be no electroweak phase transition at all || |4[][| The existence of the 
baryon asymmetry alone thus requires physics beyond what is currently known. 

The simplest extended scenario is that the baryon asymmetry does get generated at 
the electroweak phase transition, but that the Higgs sector of the electroweak theory 
differs from that in the Standard Model. In particular, it is natural to study the 
electroweak phase transition in the MSSM [0— (BJ. It has recently become clear that 
the electroweak phase transition can then indeed be much stronger than in the Standard 
Model, and strong enough for baryogenesis at least for Higgs masses up to 80 GeV ||- 



16]. For the lightest stop mass lighter than the top mass, one can go even up to 



100 GeV [IT]: in the most recent analysis [|T^], the allowed window was estimated at 



rriff rsj 75 . . . 105 GeV, mx ~ 100 . . . 160 GeV. In this regime, the transition could even 



proceed in two stages ||17|| , via an exotic intermediate colour breaking minimum. This 
Higgs and stop mass window is interesting from an experimental point of view, as well, 
as the whole window will be covered soon at LEP and the Tevatron |D| . 

Apart from the strength of the electroweak phase transition, which merely concerns 
the question whether any asymmetry generated is preserved afterwards, some dynam- 
ical aspects of the transition have also been considered recently |19||-||24|. While the 



uncertainties in non-equilibrium studies are much larger than in the equilibrium con- 
siderations, there are nevertheless indications that the extra sources of CP-violation 
available in the scalar sector of the MSSM could considerably increase the baryon num- 
ber produced, with respect to the Standard Model case |2(J . We will not consider 
the non-equilibrium problems any further in this paper, and only study whether the 
transition is strong enough so that any asymmetry possibly generated can be preserved 
afterwards. 

Even the strength of the electroweak phase transition in the regime considered is 
subject to large uncertainties. The first indication in this direction is that the 2-loop 
corrections to the Higgs field effective potential are large and strengthen the transition 
considerably flQ |. A further sign is that the gauge parameter and, in particular, the 



renormalization scale dependence of the 2-loop potential, which are formally of the 
3- loop order, are numerically quite significant fl7|l . Moreover, the experience from 



the Standard Model 0, M is that there may be large non-perturbative effects in some 



3 Unless there are relatively strong magnetic fields present at the time of the electroweak phase 
transition Q 
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regimes. Hence one would like to study the strength of the phase transition non- 
perturbatively. 

The purpose of this paper is to study the MSSM electroweak phase transition with 
lattice Monte Carlo simulations, in the regime where the transition is strong enough 
for baryon number generation. Furthermore, the results are extrapolated to the infinite 
volume and continuum limits^. Since the MSSM at finite temperature is a multiscale 
system with widely different scales from ~ nT to ~ Qw^i and since there are chi- 
ral fermions, the only way to do the simulations in practice in an effective 3d theory 
approach^. This approach consists of a perturbative dimensional reduction into a 3d 
theory with considerably less degrees of freedom than in the original theory |27|— [|30 



and of lattice simulations in the effective theory. The analytical dimensional reduc- 
tion step has been performed for the MSSM in [12|, [13|, [14], [17j. Lattice simulations in 
dimensionally reduced 3d theories have been previously used to determine the prop- 
erties of the electroweak phase transition in the Standard Model in great detail both 
at sin 2 ^ = 0, 1,0-0 and sin 2 6 W = 0.23 as well as to study QCD with 
N c = 2, 3 in the high-temperature plasma phase [39]— [41], and to study Abelian scalar 
electrodynamics at high temperatures K2], [EJ. 

The plan of this paper is the following. In Sec. |2| we formulate the problem in some 
more detail. The expressions used for the parameters of the 3d theory in terms of 
the 4d physical parameters and the temperature are given in Sec. |[ In Sec. f| we 
discuss the observables studied and the conversion of 3d results for these observables 
to 4d physical units. In Sec. | we discretize the theory, and in Sec. | we describe some 
special techniques needed for the Monte Carlo simulations. The numerical results and 
their continuum extrapolations are in Sec. [7], and in Sec. ^] we discuss the results and 
compare them with perturbation theory. The conclusions are in Sec. p]. 



2 Formulation of the problem 

In an interesting part of the parameter space, the effective 3d Lagrangian describing 
the electroweak phase transition in the MSSM is an SU(3)xSU(2) gauge theory with 
two scalar fields [|T^, [T7J : 

+ {piUyiBlU) + m 2 U3 U^U + X U3 {U^U) 2 + l3 H^HU^U. (2.1) 

Here Df = di-ig W3 t a A* and Df = di-ig S3 T A Cf are the SU(2) and SU(3) covariant 
derivatives (t a = a a /2 where a a are the Pauli matrices), gw 3 and gs 3 are the corre- 
sponding gauge couplings, H is the combination of the Higgs doublets which is "light" 



Some of the results were reported already in 25 . 
5 Finite temperature 4d lattice simulations have been performed for the SU(2)+Higgs model 
but they are extremely demanding. 
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at the phase transition point, and U is the right-handed stop field. The U(l) subgroup 
of the Standard Model induces only small perturbative contributions , and can be 
neglected here. 

The complexity of the original 4d Lagrangian is hidden in Eq. Q2.1p in the expressions 
of the parameters of the 3d theory. There are six dimensionless parameters and one 
scale. The scale can be chosen to be g§ 3 = g%T, and numerically g| 3 ~ T, see Eq. ( |3.13 ). 
Then the remaining dimensionless parameters of the theory are 



;/3 



9w3 73 = A^3 A_ 

9s3 9s3 9S3 9S3 

i — j yu — — — 

9S3 9S3 



^fl3W3j „, '"V3\yS3J (o n\ 

yu — 1 — , yu — — i — • 



Here m 2 H3 (fi), m^ 3 (/i) are the renormalized mass parameters in the MS scheme. A 
dimensional reduction computation leading to actual expressions for these parameters 
has been carried out in |T7| for a particularly simple case. Let us stress here that the 
reduction is a purely perturbative computation and is free of infrared problems. The 
relative error has been estimated in |T2], [T7|], and should be < 10%. 



The purpose of this paper is to study the theory in Eq. ( |2.1|) non-perturbatively. 
Now, the 3d non-perturbative study is completely factorized from the dimensional 
reduction step. Thus one could just study the properties of the phase diagram as a 
function of the parameters in Eq. (|2.2|) . However, a six- dimensional parameter space 
is quite large, and lattice simulations are not well suited to determining parametric 
dependences. In order to study a physically relevant region of the parameter space, 
one should hence employ some knowledge of what the reduction step tells about the 
values of the 3d parameters. At the same time, it does not seem reasonable to employ 
the full very complicated expressions. That would merely make it more difficult to 
differentiate between what is a purely perturbative effect in the dimensional reduction 
formulas and what a non-perturbative 3d effect. 

A reasonable compromise seems to be that one employs some very simple reduction 
formulas, derived in a particular special case. These are used to write the parameters 
in Eq. ( |2.2| ) in terms of fewer parameters, such as tan/5, rhu, T. The actual expressions 
used in this paper are given in Sec. || Then one studies the system non-perturbatively, 
and compares with 3d perturbation theory, employing the same 3d parameters. To 
be more precise, we compare with 2-loop 3d perturbation theory in the Landau gauge 
£ = and for the MS scale parameter ft = T, values which have been used in ||18|| , 
as well. This allows to find out whether there are any non-perturbative effects in the 
system. Once this has been done, one can go back to a more complicated situation and 
study it perturbatively, adding to the perturbative results the non-perturbative effects 
found here. Let us stress that as the reduction step is purely perturbative, the non- 
perturbative effects found with the 3d approach apply also to the effective potential 
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Figure 1: The upper bound on the 3d scalar self-coupling xh = ^^3/^53, for a fixed 
r = 0.385, xu = 0.159 (see Sec. ||), and jjh tuned such that we are at the transition 
point. This figure has been obtained with the 2- loop effective potential, requiring 
Vh/T>1. In the limit z — ► 0, it is known from |32j that the upper bound is Xh 



{^H3/ ' 9wz) r ~ O-O^ x 0.385 ~ 0.015, independent of yu. The simulation points (for the 
symmetric phase — > broken H transition) are marked with the filled circles. Since we 
have Xh = 0.0787 in the simulations (see Sec. |3|), we are always in the regime where 
the transition is strong enough for baryogenesis according to perturbation theory. 



computed in 4d [10], [16], |18|, so that the 4d potential can be used for the final studies, 
as well. 

Finally, let us note that the philosophy here has to be slightly different from that in 
the Standard Model [[52J. In the Standard Model case, there are only two dimensionless 
parameters (y = m 2 H3 /g^ 3 ,x = A^/g^), and as one of them determines the location 
of the phase transition, the properties of the transition depend only on one 3d parameter 
(x). Thus one can parameterize the non-perturbative properties of a large class of 4d 



theories in a universal way p8| . In the present case, the properties of the transition 
depend in an essential way on at least three dimensionless parameters which may vary: 
xh with the Higgs mass (or tan/3), yu with the stop mass, and z with the squark 
mixing parameters. Thus 3d lattice simulations have to be made in an essentially 
larger parameter space than for the Standard Model, and the results are in a sense 
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less universal. In particular, the constraint x < 0.03. . .0.04 [[32] for a transition to be 
strong enough for baryogenesis would be replaced by a function xh < f(yu,z). To 
demonstrate the pattern in this larger parameter space, the 2-loop perturbative results 
for f(yu, z) are shown in Fig. [I], together with the main simulation points. 

3 The parameters of the 3d theory 



To motivate the parametrization of Eq. (|2.2| ) to be given in Eqs. ( |3.13| ), consider the 
simplest case available: a large rriQ ~1 TeV, vanishing squark mixing parameters, and a 
heavy CP-odd Higgs particle {vtla > 300 GeV). Then even the formulas in Appendix A 



of |l7j can be simplified as the Q-squarks decouple: the terms involving ijiq (as well 



as mo) can be left out. The only places where one has to be somewhat careful are the 
scalar self-coupling \h where itiq re-enters in the logarithm due to zero temperature 
renormalization effects; the thermal screening terms proportional to T 2 where some 
of the contributions are to be left out due to a decoupled Q-field; and the number 
of effective scalar degrees of freedom in the gauge couplings. The results remaining 
are parameterized by tan/3 and mfj. In terms of these and m t = 170 GeV, the zero 
temperature squark masses are taken to be 

m L = -™v + m t' m h = m Q + m t + \ m z cos2/3 ; cos2/3 = j + ^ 2 ^ - ( 3 - 1 ) 

Then the Higgs mass (denoted by rhu to remind that only the leading 1-loop relation 
is employed) is 

m\ = m\ cos 2 2/? + In ( 3 .2) 

In the derivation of the finite temperature effective theory, we will take into account 
the integration over non-zero Matsubara modes, as well as the integration over the zero 
Matsubara modes of the zero components of the gauge fields A , C . With the effective 
finite temperature mass spectrum included (two SU(3) scalar degrees of freedom, U, 
D, and one SU(2) scalar degree of freedom, the Higgs doublet), the thermal Debye 
screening masses are 

11 7 
m A = ~^9wT 2 , m Co = ^9 2 S T 2 . (3.3) 

The gauge couplings appearing are approximated as 

<7|~1.1, ^-0.42, (3.4) 

according to the arguments in [|17| . These correspond to the zero temperature couplings 
ots{mz) ~ 0.12, gw( m z) ~ 2/3. The effects of the U(l) gauge coupling are small and 
will hence be approximated by g' — 1/3. 
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From these approximations and the formulas in 17], it then follows that 

9 2 W T 



2 

9\V3 



9S3 
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73 
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1 



247rmA 
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lGir 
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1 
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h 2 sin 2 (3T, 



(3.5) 

(3.6) 

(3.7) 

(3.8) 
(3.9) 



~ 2 

m 



^ ( ^(9 2 w + 9' 2 ) cos 2 2P + + l^' 2 + i/^ 2 sin 2 /3 ) T 2 



-g w Tm Ao 



+ 



H3 



+ 



16tt 

16^2 ( ^9W3 + ^H39W3 - 12X H3 ~ 3 7s + 8 9S3ls) ln ~- 
2 , ' 1.2 . ' ' 

' 8 9s3 + Y Xu39 S3 ~ 16X us ~ 2 7s + 3^73) In ■ 1 " 



+ ( n9 S + a h t Sin P ) T - ^9s Tm C 



167T 2 

,2 • 2o _ „2 ™2 //r> 2 



(3.10) 



(3.11) 



Here /i 2 sin /5 = g^m 2 / (2171^) . In [I7[], it was argued that as a first estimate, 

Ai/3 ~ At/3 ~ 7T, (3.12) 
and we will use this assumption here. The precise determination of Ah3, ^U3 requires 



a 2-loop dimensional reduction computation, and has recently been carried out in ||44 
Using eqs. ( |3.5| )-( |3TT2] ), we finally get a simple parametrization for Eq. 



9S3 

r 

z 

X H 

Vh 

yu 



1.085T, 

0.385, 

0.893, 

0.159, 

yi(m H ) - y 2 (rh H ) 



AOOGeVN 2 



U.517- 0.819 



V T 

2 



(3.13) 



Here the tan/3 (or rhjj) -dependent functions are given in Table |I[ In what follows, we 
study the 3d theory in Eq. ( |2.1| ), parameterized by T,m,H,rhu through Eqs. ( |3.13| ). 
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fn-nft 
IdbLLfJ 




f\ m H) 




V2 \Jlt>H) 


3 


94.5 


0.0787 


0.548 


0.379 


5 


103.4 


0.0917 


0.554 


0.454 


7 


106.1 


0.0960 


0.556 


0.478 


9 


107.3 


0.0978 


0.557 


0.489 


12 


108.1 


0.0991 


0.558 


0.496 


20 


108.8 


0.1002 


0.558 


0.502 


30 


109.0 


0.1005 


0.558 


0.504 



Table 1: The functions appearing in Eq. ( |3.13| ). 



To display the general phase structure of the theory, consider the 2-loop effective 
potential JDJ- For definiteness, we consider the Landau-gauge and the MS scale pa- 
rameter p, = T. The phase structure following from the 2-loop potential is shown in 
Fig. |2|. The general pattern is that when the 3d parameters are as in Eq. ( |3.13|) , the 
system has a first order transition at T c ~ 100 GeV for rhy < 60 — 70 GeV. This tran- 
sition is rather strong even though fhu is large, due to the stop loops. As rhu becomes 
larger {mi smaller), the transition gets even stronger, and then at some point one may 
get a two-stage transition |T7[ . The existence of a two-stage transition depends on the 
parameters of the theory, and for large squark mixing parameters the parametrization 
in Eq. ( |3.13|) would change so that the two-stage region is not reached |L8|. It will 
be seen below that the qualitative behaviour in Fig. |2| is reproduced by Monte Carlo 
simulations. 



4 Observables in 3d and 4d 

Before going to simulations, let us discuss how the physical observables to be measured 
in 3d are related to the corresponding 4d observables. 

Consider first the Higgs field vacuum expectation value vh- This is the object by 
which one usually characterizes whether the phase transition is strong enough for baryo- 
genesis |], |J, the requirement being vh/T> 1. As such vh is, however, a gauge de- 
pendent quantity. If one computes it in the Landau gauge (vjj), as is usually done, 
then in terms of gauge-invariant operators the same expression would be non-local. On 
the other hand, there is a simple local gauge-invariant quantity closely related to %, 
namely WH ~ v 2 H /2. The problem with H^H is that being a composite operator, it 
gets renormalized beyond tree-level. One can consider two recipes for defining a unique 
scale independent object: either one considers the discontinuity of H^H between the 
broken and symmetric phases, or one considers, say, the MS-regularized WH at the 
natural scale g| 3 . The first option differs qualitatively from the Landau-gauge in 
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Figure 2: The perturbative phase structure using the parametrization in Eq. ( |3.13|) , 
together with the 2-loop Landau gauge vev v#/T in the broken phase. 

that given T c , it is sensitive to the value of WH in the symmetric phase and is thus an 
inherently non-perturbative quantity, contrary to vjfr. Moreover, it can only be deter- 
mined at the phase transition point (or the metastability region) where a discontinuity 
exists, unlike v^. We thus choose the latter option which gives, in principle, a purely 
perturbative broken phase quantity for given T. Hence we define 

vh_ 
T 

vu_ 
T 

Note that due to a trivial rescaling with T, the dimension of H^H, U^U is GeV in 3d. 

Other observables needed in the study of the phase transition are the latent heat, 
the interface tension, and the different correlation lengths. These enter, for instance, 
the estimates for the nucleation and reheating temperatures (see, e.g., |[7|]); which are 
needed to decide whether vh/T should be taken at T c or some other temperature. 

Consider first the latent heat. It is determined by 

L = T c -^Ap(T) T=Tc , (4.2) 



2^r< 






\ 9S3 






T ( 


\ 9S3 



1/2 



1/2 



(4.1) 



8 



where Ap(T) = p s (T) — pb(T) and p s ,Pb are the pressures in the symmetric and broken 
phases, respectively. At the phase transition point, the partition function obeys 

| 4 = d exp(^Ap) ^= d exp(-V<4A6 3 ), (4.3) 

where €3 is the dimensionless vacuum energy density of the 3d theory in Eq. ( 27T ) . 
As such, €3 is a divergent quantity, but its jump across a phase transition is finite, 
and the derivatives of the jumps with respect to different parameters are related to 
jumps of different operators. In particular, as only the parameters yn,Vu depend on 
the temperature according to the parametrization in Eq. ( |3.13|) , one gets 

^4 = -j^—Ae^Z^XH^Xu^HM 



- c 



9k (dy H /gtg\ dy v 




= f (f a(^Wa ( ^)] , (4.4) 

where A (...> = -A (. . .) = (.. .> 6 - (. . .>.. 

Consider then the interface tension. The interface tension is the extra free energy of 
a phase boundary separating two coexisting phases. Due to this extra free energy, the 
probability P m i n of a configuration where there is such an interface, is smaller than the 
probability P max for a pure phase. In 4d and 3d units, this is expressed as 

P ■ 1 

j^- = exp(--o-A) = exp(-o- 3 A^ 3 ). (4.5) 

-'max 

It follows that on a lattice with periodic boundary conditions and the cross-sectional 
area A = (N x a) 2 , 



mm 



where it was taken into account that the boundary conditions force there to exist at 
least two interfaces. According to Eq. Q4.5Q , the physical interface tension is then 

£ = f?*3 = (1.085) V 3 . (4.7) 

c c 

Finally, let us discuss correlation lengths. The lowest-dimensional gauge invariant 
continuum operators, used for the mass measurements, are 

1. 



S H = WH, G H = -/•;"/',, 

V l H = ImHWfH, 



4 
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Note, in particular, that while there are perturbatively even massless gauge excitations 
in the phase with a nonzero vu, these are confined by an unbroken SU(2) subgroup and 
only the SU(3) singlet excitation corresponding to Vjj is physical [fi5]| . The real parts 
of WDfH, WD?U could be used, as well, but they are expected to couple to the same 
excitations as Sh, Sjj. There can also be non- vanishing cross correlations between the 
operators, e.g. between Sh, Sjj, and in such a case the true eigenstates can be obtained 
by diagonalizing the correlation matrix. The correlation lengths are measured in units 
of g| 3 , and can then be trivially converted to units of T. 



5 The lattice action 

We now discretize the theory in Eq. ( |2.1| ) with standard methods. The scalar fields 
are rescaled into a dimensionless form by H^H — > H^H(3H/{2a) = WHg 2 S3) U^U 
WU/3u/(2a) = WUgg 3 . The lattice Lagrangian is then 

Matt 



i<j 



1 



-ReTri^f 



1 



-Tr P w 

2 13 . 



- fa E Re # f (x)^T + i) + $iPH + P? (H^H) 2 

i 

- (3uJ2 Retr f (x)£/f (x)C/(x + i) + $Wu + (3% (fffjf 

i 

+ PJH^HWU. (5.1) 

Here Pfj,Pff are the SU(3) and SU(2) plaquettes, respectively, and Uf(x), U^(x) are 
the corresponding link matrices. The SU(2) Higgs field is conveniently expressed with 
a 2 x 2 matrix parametrization: 



H 



H* H x 
-H* H 2 



ho&0 + i E haVa, 



(5.2) 



0=1 



where <Jq = 1 and u a >o are the Pauli matrices. When we write the Lagrangian (15 
in terms of the representation (|5.2p , the ff-terms are substituted through WMH — > 
|Tr H^MH. This is the form of the action we actually use in the simulations. 

The lattice parameters appearing in Eq. ( |5.1|) can be expressed in terms of the lattice 
spacing a and the continuum parameters in Eq. ( j2.2| ). Let us denote 

0s = -4-- (5.3) 
a 9ss 

Then it follows from the discretization procedure and from the relation between the 
MS and lattice regularization schemes in 3d |I6| , |47|j that 

2p s „ 216 
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36 216 f 

+ 



Ps PI V 



yu - - + 8x v + 2z 



^Ps 
4tt 6 



1 



64 



8 + ^ x u ~ 16a& + 3rz - 2z 2 ) (In p s + 0.08849) 



16vr 2 

+ 19.633 + 12.362xc/ + 1.7384rz 



(5.4) 



Thus all the lattice couplings are determined, once the continuum parameters and 
Ps have been fixed. The relations in Eq. ( |5.4p become exact in the continuum limit. 
Improvement formulas at finite lattice spacing which should remove the 0(a) effects 
from most of the quantities have been derived in (see also |49[|). 

To measure the observables in Eq. ( |4.1| ), one needs the relations of the MS and lattice 
regularization schemes also for H^H, U'U. It follows f47| that 



53) 



9ss 



WH 



12vr 



Ps 



16tt 



; r(lnp s + 0.66796), 



/WU(g 



2 ) 

S3) 



9S3 



U ] U) ~ ^Ps - t^^Ps + 0.66796). 



(5.5) 



Note that the discontinuities of H^H, U^U are finite. 

With these relations fixed, we are ready to go to simulations. Extrapolations to the 
infinite volume and continuum limits will then allow to determine non-perturbatively 
the properties of the MS continuum theory in Eq. fl2.1|) . 



6 The Monte Carlo update algorithm 

There are three basic reasons which make the lattice simulations of the theory in 
Eq. ( |5.1| ) quite a demanding numerical problem. 
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1. For the simulations to be reliable and allow an extrapolation to the infinite volume 
and continuum limits, the lattice spacing a and the smallest linear extension L = Na 
of the lattice must satisfy 

a < £min < £max < Nd, (6.1) 

where £ min and £ max are the smallest and the largest physical correlation lengths in 
the system. Thus a multiscale system where £ m i n <C £ ma x, requires very large lattice 
sizes N. The present system does have many different excitations and scales, part 
of them proportional to g§ 3 and part to g^r 3 <C gf 3 . We will measure the different 



correlation lengths in Sec. |7]5|. It turns out that Eq. (|6.1j) can actually be well enough 
satisfied, due to the fact that the transition is quite strong. In a very weak (or second 
order) transition, the determination of (non-universal) observables would require much 
larger lattices, since some of the correlation lengths are very large. The lattices used 
are shown in Table [| 

2. While a strong transition makes it easier to satisfy Eq. ( |6.1|) , there is at the same 
time a serious new problem. Indeed, the transition can become so strong that during 
the Monte Carlo simulation the system does not want to tunnel from one metastable 
minimum to the other, especially for the large volumes needed in order to satisfy 
Eq. ( |6.lD . On the other hand, one needs to probe both phases simultaneously with 
sufficient statistics in order to reliably determine the relative weights of the phases 
(important for determining the transition temperature) and the suppression of the 
mixed phase (necessary for interface tension measurements). To allow for sufficient 
tunnelings, one has to use multicanonical simulation algorithms. 

3. A special feature of the Lagrangian ( |5.1| ) is that in the large fhu -region both the 
H and U fields play a significant role in the transition, and both fields can become 
"broken" , though not at the same time. The coupled dynamics of the H and U fields 
makes the optimization of the update algorithm a delicate issue: it is only too easy 
to select an update move which evolves the configurations through the phase space 
extremely slowly. This is especially relevant in multicanonical simulations, where, as 
we shall see, the choice of the multicanonical order parameter becomes critical. 

In the following subsections, we discuss in some detail the methods employed to meet 
these (partly exceptional) requirements. 



6.1 The overrelaxation update 

As usual in simulations of a system which undergoes a phase transition, the overrelax- 
ation update is much more efficient in evolving the fields than the diffusive Metropolis 
or heat bath updates. Intuitively this is easy to understand: the overrelaxation update 
propagates information through the system in wave motion (dist. oc t), whereas heat 
bath and Metropolis obey the diffusion equation (dist. oc t 1 ^ 2 ). However, in order to 
ensure ergodicity one has to mix heat bath -type updates with overrelaxation. 
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We use a compound update step which consists of 4-6 overrelaxation sweeps through 
the lattice followed by one heat bath/Metropolis update sweep. In one sweep we first 
update all of the gauge fields, followed by the updates of the Higgs fields. 

The gauge field update. Compared with the Higgs fields, both SU(2) and SU(3) 
gauge fields are relatively 'inert' with respect to the transition, i.e., their natural modes 
evolve much faster than the Higgs modes (the slow gauge modes arise only through the 
coupling to the Higgs fields). Thus, the gauge field update algorithms are not as critical 
as the ones for the Higgs fields. We use a gauge field update not qualitatively different 
from the standard SU(2) and SU(3) pure gauge updates, in spite of the hopping terms 
of the form Tr $' (x)C/j(x)$(x+z) in the action (note, however, the modifications due to 
the multicanonical update, Sec. |6.2|) . We use the conventional reflection overrelaxation 
and Kennedy-Pendleton heat bath |5U| methods; for SU(3) the updates are done on 
the SU(2) subgroups of the SU(3) link matrices. 

The overrelaxation of the Higgs fields. Efficient overrelaxation update algo- 
rithms for the Higgs fields H and U are essential in order to minimize the autocorrela- 
tion times. From the lattice Lagrangian ( p.lp we can observe that the local action for 
the Higgs fields H(x) and U(~x) can be written in the following generic form: 

V[(f>(x)] = -F o (x)0 a (x) + C 2 (x)i? 2 (x) + C 4 i? 4 (x), (6.2) 

where (f> is either U or H, R 2 = <f) a <p a , and the index a is understood to go through the 
real and imaginary parts of the components of the complex vectors separately: thus, 
for H, a = 1 ... 4, and for U, a = 1 ... 6. F(x) is the sum of the nearest-neighbor 
'force' terms at site x: 

F{x) = /3 H Yl [^(x-2)0(x-O + ^(x)0(x + z)]. (6.3) 

1=1,2,3 

This form seems to suggest separate update steps for the radial and SU(iV)-components 
of the Higgs fields. However, as already noticed in the simulations of SU(2)+Higgs 
systems in 3 and 4 dimensions [^B], [32|, it is much more efficient to perform an update 
which simultaneously modifies the radial component and the direction of the Higgs 
fields. 

In our Higgs field update we generalize the Cartesian overrelaxation , presented in 
Ref. we update the Higgs variables in the plane defined by 4- or 6-dimensional 
vectors 4> a and F a , using the components of (j) a parallel and perpendicular to F a : 

X = f a <j) a} Y a = <f) a -Xf a , (6.4) 

where f a = F a /F and F = ^F a F a . In terms of X and Y a , Eq. ( |6.2|) becomes 

V(X, Y) = -XF + C 2 (X 2 + Y 2 ) + 2C 4 X 2 Y 2 + C 4 (X 4 + Y 4 ) . (6.5) 
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The overrelaxation in Y is simply the reflection Y a — > — Y^, or a — > —<f) a + 2Xf a . For 
SU(2), this is exactly equivalent to the conventional reflection overrelaxation procedure. 
For the X-component, the polynomial form of V(X) provides a way to perform an 
efficient approximate overrelaxation: we find the solution to the equation V(X') = 
V(X) and accept X' with the probability 

p m=mi n(p ,l), P. = ^P. (6-6) 

Since V(X) is a fourth order polynomial, solving the equation V(X') = V(X) reduces 
to finding the zeros of a third order polynomial (we already know one zero X' = X, 
which can be factored out). The parameters of V(X) are such that there always is 
only one other real root, and it is straightforward to write a closed expression for X'. 
The update is an almost perfect overrelaxation: in our simulations the acceptance rate 
varies between 99.4% - 99.9% for both U and H, depending on the fis used. The 
acceptance is high enough so that the "diffusive" update dynamics inherent in the 
Metropolis accept/reject step does not play any role, and the evolution of the field 
configurations is almost deterministic. 

The essential part of the Cartesian overrelaxation is the X-component update; the 
update of Y has only a small effect on the evolution of the fields. Intuitively, the X- 
mode update achieves its efficiency by suitably balancing the entropy and the action: 
in a single update move, it interpolates between states where (i) |<3>| is large and the 
direction is relatively parallel wrt. its neighbours and (ii) |<&| is small and the direction 
more 'randomized'. Update steps acting separately on the length and the direction of 
the Higgs fields do not achieve this kind of balancing, and hence the magnitude of the 
change in a single update can be much smaller. 



6.2 The multicanonical update 

The first order phase transitions are relatively strong in the whole parameter range 
studied in this paper. Thus, in standard simulations using the canonical ensemble, the 
tunnelling rate from one phase to another becomes very small at all appreciable lattice 
volumes. This probabilistic suppression is due to the existence of the phase interfaces 
in the mixed phase, and it is proportional to the interface tension times the area of the 
interfaces (see, for example, Fig. ^). 

To enhance the probability of the mixed states we can modify the lattice action with 
the global multicanonical weight function W: 

Sue = E Aatt(x) - W(R) , R = Y1 r(x) , (6.7) 

X X 

where r(x) is a suitable local order parameter sensitive to the transition. The canonical 
expectation value of an operator O can be calculated by reweighting the individual 
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multicanonical measurements Ok with the weight function: 

(O) = J2 O k e~ w ^/J2 e_WW , (6-8) 

k k 

where the sums go over all measurements of O and R. In this work we use a continuous 
piecewise linear parametrization for W{R) (see Eq. ( |6.1(J| )). 

Selection of the multicanonical variable. The structure of the phase diagram 
of the theory in Eq. ( |5.1|) is complicated enough so that the choice of a suitable multi- 
canonical order parameter r(x) is, in general, not obvious. In the small rhu -region the 
transition is driven by the SU(2) Higgs field H, and a good choice is r = wH. This 
order parameter has been widely used in standard SU(2)+Higgs simulations |3"B 



However, the situation is very different when rhu becomes so large that one is near 
the triple point in the phase diagram (see Fig. ^)): the SU(3) Higgs field U becomes 
important, and WH is not necessarily an optimal order parameter any more. Indeed, 
it turns out that the selection of an optimized multicanonical variable is essential for 
the performance of the algorithm. In the following we discuss the different choices used 
for the variable. 

In Fig. |3] we show the joint probability distribution of WH and WU near the triple 
point, i.e., the coexistence point of the symmetric, broken H, and broken U phases. The 
distribution is from an rhu = 67.05 GeV, T = 84.3 GeV, volume 12 3 , fig = 12 lattice. 
The appearance of three peaks in the probability distribution is clearly visible. The 
peak corresponding to the broken H phase is strongly separated from the symmetric 
phase and broken U phase peaks, whereas there is only a mild suppression between 
the broken U and symmetric phase peaks. This is a clear signal of the strong first 
order transition between the broken H and symmetric phases, and relatively much 
weaker transition between the broken U and symmetric phases. The broken phases are 
connected only through the symmetric phase. 

From this figure one can already see that WH does not distinguish the symmetric 
phase and the broken U phase. Indeed, on the left panel of Fig. |] we plot the one- 
dimensional probability distribution p(WH) from the data shown in Fig. || The 
symmetric phase and the broken U phase fall on the same peak in the distribution. Thus 
one should use some other variable which can distinguish the whole phase structure, 
in order to enhance also tunnellings from the broken U phase to the symmetric phase. 

Motivated by Fig. ^ an obvious choice would be to use a weight function of two 
variables, W(WH,WU). While this would certainly be possible in principle, the use 
of a weight function with a higher than one-dimensional argument is quite cumbersome 
in practice and we did not attempt to do it here. Instead, we shall consider the one- 
dimensional weight function variables (H^H — WU) and (Vg — Vu), where Vh and 
Vu are the Higgs field "hopping terms" , where the length of the Higgs fields has been 



15 




Figure 3: The emergence of the triple point: the joint probability distribution of the 
SU(2) and SU(3) Higgs field lengths squared, H*H and W&, at m v = 67.05 GeV, 
T = 84.3 GeV, on a 12 3 , f3 s = 12 lattice. Here (O) = ^ x O{x)/V. The three peaks 
correspond, from left to right, to broken H, symmetric, and broken U phases. The 
relative strength of the transitions is evident from the suppression of the probability 
density between the peaks. When the volume is increased, the suppression between 
the peaks grows and the peaks become sharper. 

divided out: 

V h = ^7E Re/it ( x ) f/ r( x ) /l ( x + 2 ) H = R H xh, \h\ = l, 

1 x > 1 (6.9) 

Vu = — ^Rew^xjf/f (x)u(x + i) U = RuXu, \u\ = 1. 

3^ x ,i 

In contrast to the distribution p(WH), the distribution p(V# — Vu) on the right panel 
of Fig. H clearly separates the three phases. We thus expect (Vh — Vu) to be a much 
better multicanonical order parameter than (the volume average of) WH. Indeed, in 
Fig. [5] we compare the performances of three multicanonical algorithms with weight 
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Figure 4: One-dimensional probability distributions for the system shown in Fig. |3|. 
Left: the distribution of WH = J2 X H^H(^)/V. The left peak corresponds both to 
the symmetric phase and the broken U phase, the right peak to the broken H phase. 
Right: the probability distribution of (Vh — Vu), where Vh and Vu are the SU(2) and 
SU(3) hopping terms without the Higgs field radius, Eq. fl6.9|). The three phases are 
clearly separated: from left to right, the peaks correspond to the broken U, symmetric, 
and broken H phases. 



functions W(R) and variables R = H^H (top panel), R = WH — WU (middle panel), 
and R = Vh — Vu (bottom panel). The time histories are from simulations of an 
rhu = 68 GeV, (3s = 12, 8 3 -system at T = 83 GeV. At these parameter values the 
system has a strong first order transition between the broken H and broken U phases. 
It is evident that the last multicanonical algorithm is much more efficient than the 
two first ones in driving the system through the transition. It should be noted that 
the volume 8 3 used here is truly microscopic and useless for a quantitative analysis; in 
any reasonable volume we had trouble to make the first two algorithms to tunnel even 
once. 

The multicanonical variable (Vh — Vu) depends on all of the fields on the lattice, 
and it has to be evaluated after each stage of the update. Why is it better than 
the formally simpler variable (H^H — WU)? One possible explanation is the fact 
that the probability distributions p H (H^H) and pu(WU) are, by themselves, very 
asymmetric: they have a sharp symmetric phase peak and a broad broken phase peak 
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Figure 5: Sections of the Monte Carlo time histories of the observable Vh from simula- 
tions where the multicanonical order parameter is WH (top), (WH — WU) (middle) 
and Vh — Vv (bottom). 



(see Figs. [|, [] and 0, for example). Remembering that in the broken U -phase, H is 
always in the symmetric peak and vice versa, and recognizing that around any single 
peak the probability distribution p H -u(WH — UW) is in essence a convolution of 
Ph and pu, we see that the distribution pH-ui.%) has broadened peaks for both the 
broken U and broken H phases. Thus, the power of the multicanonical weight function 
W(WH - WU) is 'softened'. 

This can be contrasted to the probability distributions p(Vh) and p(Vu) (Fig. |]). 
Now both the symmetric and the broken peaks are of comparable width, and the 
convolution does not soften the structure of the peaks too much. 

In the simulations in this paper we use two different multicanonical variables: H^H 
in the small rhu -region, and (Vh — Vj) when fhu is large. Since we are using a parallel 
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supercomputer, it is not practical to keep track of the value of the global multicanonical 
variable when the individual Higgs and gauge field variables are updated. We calculate 
the value of the multicanonical weight variable only after each global even/odd -site 
update sweep, and perform an accept/reject step for the whole update. Nevertheless, 
the change in the multicanonical variable remains small enough so that the acceptance 
rate is better than 90-96%, depending on volume and (3s- 

Recursive calculation of the weight function. When the system has a first 
order phase transition, the goal is to choose the weight function W(R) such that the 
resulting probability distribution puc(R) is approximately constant in the interval 
Ri < R < R2, where R\ and R 2 denote the pure phase peak locations. This is one 
of the main difficulties of the multicanonical method: a priori, the weight function 
is not known, and the optimal weight function is W(R) = — lnp can (_R), where p can 
is the canonical probability distribution of the observable R. This is one of the very 
quantities we attempt to determine with Monte Carlo simulations. Thus one has to 
use some sort of an iterative procedure for determining W(R). 

A precise determination of W(R) is needed especially in the regime of large fhu, 
where the first order transition between the broken H and U phases becomes extremely 
strong. If we allow that the multicanonical probability distribution is 'ideal' up to a 
factor of, say, 1.5, then the multicanonical weight function must be determined with an 
absolute accuracy In 1.5 pa 0.5. The largest variation of W in this study is ~ 100. (That 
is, we have to boost the probability of the mixed state with respect to the pure phases 
by a factor pa exp 100.) Thus, the weight function has to be determined to an overall 
relative accuracy of 0.5%. We determine the weight function with an automatized 
recursive process.^] 

We parameterize W with a piecewise linear continuous function: 

W(R) = Wi + (w i+1 - Wi ) J* - % , Ri < R < R i+1 . (6.10) 

The fcth estimate of the weight function is W k , and W l is set to the initial estimate (it 
can also be a constant function). The iterative process we use here to improve on the 
estimate wf is based on the relative (canonical) probabilities Pi that the system is in 
bin (i — 1) or i: then the weights are chosen so that Wi — iUj_i = ln(pj_i/pj). In more 
detail, the method proceeds as follows:[] 

(i) During a (short) run of M iterations, measure 

M M 

n\ = £ ^{R m ) and h\ = £ e~ wk ^5 t (R m ), (6.11) 

m=l m=l 

6 A somewhat different recursive method for calculating W has been described in Ref. [|5l| . 
7 In this discussion we ignore the special treatment required by the boundary bins and by bins of 
unequal width. 
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where Si(R) = 1, when R4-1/2 < R < R4+1/2, otherwise 0. Thus, n\ is the number of 
'hits' in the bin number i, and h k is the measured estimate of the canonical histogram, 
(ii) After k runs, W k+l can be obtained by calculating 
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where the function g\ is a suitably chosen two-bin weight factor, characterizing the 
statistical importance of the run k. We use here gf = (nf„ 1 + n k ), if 
else g k = 0. The number n m i n ~ 10 guarantees that the bins have some minimal 
amount of statistics before they are taken into account in the calculation. The initial 
weight function W 1 can be included in Eq. ( |6.12| ) by setting h® = exp(-wj), and g® to 
a constant value proportional to the estimated 'quality' of W 1 . 

(iii) In practice, the convergence can be greatly accelerated by an overcorrection 
of W: let us calculate a modified weight function 



with a suitably chosen constant C ~ 1. Using now instead of wf in the next round 
of simulations in step (i), regions of the phase space not yet frequently visited (small 
J2 n i) are strongly favoured. This can dramatically accelerate the initial coverage of 
the whole R -range of interest, and guarantees a rough estimate of the final weight 
function only after a modest number of iterations. Naturally, the estimate of the true 
weight function (which is used in the final simulations) is still given by Eq. (16.12 ). 

(iv) The process (i)-(iii) is repeated until a good enough convergence for W k has 
been obtained. 

In our simulations, the length M of the runs in step (i) was 500-2000 iterations, 
depending on the volume, and the process (i)-(iii) was repeated ~ 100 times (total 
of 50000-150000 iterations). The whole procedure is automatized, except for the 'exit 
condition' in step (iv). 

At the beginning of these set-up simulations the changes in the weight function W 
are quite large and the system does not reach any kind of an approximate equilibrium 
distribution. When the run progresses, the amplitude of the modifications to W de- 
creases smoothly. These preliminary runs were only used for determining W, and were 
discarded for the analysis described below. 

7 Simulations and results 

Since our main interest is the phase diagram and the observables which quantify the 
strength of the transitions, most of our simulations were performed at and immediately 
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Table 2: The lattice sizes and spacings used. Multicanonical simulations are marked 
with the subscript ( m ). 



around the phase transition parameters. The lattice sizes used are listed in Table |2|. 
For each lattice listed we performed 60 000 - 200 000 compound iterations (4 x overre- 
laxation + 1 x heat bath). Some of the points shown include several separate runs at 
slightly different T, which are then combined with the multiple histogram reweighting. 

A well controlled infinite volume limit is essential in order to obtain reliable results 
for quantities like the latent heat and interface tension. Thus, we always perform 
simulations with several lattice volumes at any given lattice spacing. A cylindrical 
lattice geometry is needed especially for the interface tension, and most of the lattices 
in Table ^| are highly asymmetric. 

To extrapolate to the continuum limit, we have made simulations with two different 
lattice spacings, (3$ = 12, 20, at ifiu = 50 and 65 GeV. This only allows a linear 
extrapolation. However, it is understood analytically that the dominant corrections 
are linear |4"8"|, |4T| , and moreover, linear extrapolations work extremely well for the case 



of the Standard Model 0, pSfl . Note also that for the SU(2) coupling the /?s's used 
correspond to (3w ~ 21, 35, which are larger than the largest inverse lattice spacings 
used in |32], |33], Q . We are thus confident that the linear extrapolations provide good 
estimates of the continuum values. Moreover, as we shall see, in several observables 
the lattice spacing dependence is surprisingly small. 
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The Monte Carlo simulations were performed on a Cray T3E parallel computer at 
the Center for Scientific Computing, Finland, using 32-64 nodes. The performance of 
the code was ~ 105 Mflop/ second/node. The parallel communication parts of the code 
are based on the MILC collaboration public domain lattice QCD code p2]| . The total 
cpu-time used was about 7.5 cpu- years of a single node's capacity, corresponding to 
~ 2.5 x 10 16 floating-point operations. 

7.1 Phase diagram and critical temperatures 

The critical temperature can be determined accurately from the Monte Carlo data. As 
in the standard electroweak model, there are no known local order parameters, which 
would acquire a non-zero value only in one of the phases of the model. Instead, we use 
quantities which display a discontinuity at the transition points (when V — > oo). The 
quantities we use are WH, WU, and the hopping terms Vh and Vu, Eq. (|6.9| ). 

For each individual lattice listed in Table |2| we locate the pseudocritical temperature 
T c with several different methods (cf . Ref . |H| ) : 

(1) maximum of the susceptibility Xh\h = 

(2) maximum of the susceptibility xv H = V ((Vh — (Vh)) 2 ), 

(3) "equal weight" T- value for the distribution p(H^H), 

(4) "equal height" T- value for the distribution p(V#). 

The items above are for transitions between the symmetric phase and the broken H 
phase. For the symmetric <-> broken U transitions we use the corresponding operators 
where H — > U . For the broken U <-> broken H transition we can use all of the above 
operators and combinations thereof. The T c - values corresponding to the above criteria 
are found with the Ferrenberg-Swendsen (multi) histogram reweighting [|53j| , and the 
error analysis is performed with the jackknife method, using independent reweighting 
for each of the jackknife blocks. 

In Figs. |] and |7| we show examples of pseudocritical histograms at fhjj = 65 GeV 
and 68 GeV. At fhjj = 65 GeV, the transition is between the symmetric phase and 
the broken H phase. At rhjj = 68 GeV, there are two transitions: from the symmetric 
phase to the broken U phase (bottom panel on Fig. [7]) and from the broken U phase to 
the broken H phase (top panel). The suppression of the probability between the peaks 
characterizes the strength of the transitions: the symmetric <-> broken H transition is 
a relatively strong first order transition; the symmetric <-> broken U transition is much 
weaker; and the broken U <-» broken H transition is extremely strong. 

The infinite volume and continuum limits. For any given lattice the criteria 
(l)-(4) above yield different pseudocritical temperatures T c . However, in the thermo- 
dynamic limit V — ► oo all of the methods extrapolate very accurately to a common 
value. This is shown in Fig. || for my = 50 GeV. It should be noted that the differ- 
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Figure 6: The pseudocritical histograms at rhu = 65 GeV, (3 s = 12. Top: the equal 
weight histograms p(WH). Bottom: the equal height histograms p(Vh)- 



ent methods do not give statistically independent results, and combining the results 
together is not justified. For definiteness, we use the criterion (3), the equal weight 
of W H -histograms (or W [/-histograms, where appropriate), to determine our final 
results for T c . In strong first order transitions the equal weight criterion is very robust, 
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Figure 7: The equal weight histograms at fhu = 68 GeV, (3s = 12. Top: histograms 
of WH at the broken U <-> broken if transition. Bottom: histograms of f/W at the 
symmetric <-> broken [/ transition. 

and yields practically identical results for all suitable order parameters. 

The results of the infinite volume extrapolations for all the different parameter values 
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Figure 8: The infinite volume extrapolation of the critical temperature at mjj = 50 GeV 
for (3 S = 12 and 20. 

are shown in Fig. |]. For rhu = 50 GeV and 65 GeV, we have data with two different 
lattice spacings, and a continuum extrapolation linear in l/(3s = a 9si,/^ ls possible. 
The results, together with the perturbative results, are also shown in Table |3|. We 
discuss the comparison with perturbation theory in more detail in Sec. 

The two different lattice spacings do not offer the possibility to check for subleading 
corrections in l//3g. However, as already discussed in the beginning of this Section, the 
results from the electroweak simulations |32], [38| strongly suggest that the linear term 
dominates the extrapolation. It is also evident from Fig. || that the variations from 
(3s = 12 to the continuum limit are much less than the difference from the perturbative 
results. Thus, even the (3$ — 12 results give a fair estimate of the continuum critical 
temperatures. 

7.2 Scalar field expectation values 

The scalar field expectation values vh/T (in the broken H phase) and Vjj/T (in the 
broken U phase) are calculated from Eqs. Q4.1]) and fl5.5|) . For this we need (WH) 



results at T c by rejecting the symmetric phase and mixed phase measurements. This is 




These are calculated from the multicanonical 
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Figure 9: The phase diagram and the critical temperatures. The continuous lines are 
from 2-loop perturbation theory in the Landau gauge. 



done by imposing a suitably chosen lower cut-off for the measurements of H^H or U'U. 
For example, at rhu = 65 GeV and (3s = 12 (Fig. ||), we accept only values WH > 1.85 
(in other words, we simply calculate the center of gravity of the WH > 1.85 -part 
of the histogram). Since the mixed phase is very strongly suppressed, the V — > oo 
extrapolations are quite insensitive to the value of the cut-off. 

We calculate Vff/T and Vu/T for each of the volumes separately, and extrapolate to 
the infinite volume linearly in the inverse (smallest) cross-sectional area of the lattices. 
As in the electroweak case |32], |38| , the results from lattices of different geometries obey 
the l/A-behaviour much better than the inverse volume law. As an example, in Fig. [H] 
we plot vh/T from rhu — 50 GeV lattices. The V —>■ oo extrapolations are shown in 
Fig. [12] and in Table [3] together with the perturbative results. For rhu = 50 GeV and 
65 GeV, we again perform the continuum extrapolation linearly in l//3s = ag S3 /Q. 
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Figure 10: The infinite volume extrapolation of the H field expectation value vh/T at 
T = T c in the broken phase, for rhy = 50 GeV, fis — 12 and 20. 
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Figure 11: The infinite volume extrapolations of vh/T and vy /T at T C) 2, and vyjT at 
T Cj i, for rhy = 68 GeV. In each case the system resides in the relevant broken phase. 
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Figure 12: The scalar field expectation values in the broken phase at T c . 



7.3 Latent heat 

The latent heat can be calculated from Eq. ( |4.4| ). However, in order to reduce statistical 
noise in the measurements, it is worthwhile to transform Eq. (|4.4|) into 



™2 



™2 



T 4 T, 



7^2 

± c 



(7.1) 



Here A (. . .) = (. . .) brokon — (. . -) symm • Eq. (|7.1|) implies that we measure the symmetric 
and broken phase expectation values of the whole expression in the angular brackets, 
instead of doing it for H^H and WU separately. 

As a function of the lattice volume and spacing, the latent heat behaves very much 
like the scalar field expectation values vh and Vy. We show the results of the V — > oo 
and a — > extrapolations in Fig. 13 and in Table S 



7.4 Interface tension 

As described in Sec. |], we measure the interface tension with the histogram method. 
Due to the extra free energy of the phase interfaces in the mixed phase, the probability 
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Figure 13: The infinite volume values for the latent heat. 



of the mixed phase is suppressed by a factor oc exp(— / int /T) = exp(— aA/T). This 
is seen as a 'valley' between the peaks corresponding to the pure phases in the order 
parameter histograms; see, for example, Figs. || and [F[ 

For the interface tension measurements it is advantageous to use lattices of cylindrical 
geometry, L z 3> L x = L y . Because of the periodic boundary conditions there are at 
least two interfaces which span the lattice, and because of the cylindrical geometry, 
they tend to form parallel to the (x, y) -plane. L z should be then long enough that 
the two interfaces do not interact appreciably: this is seen as a flat minimum in the 
histograms. 

We used the histograms p(Vn) (Eq. (|6.9| )) to extract the interface tension. These his- 
tograms have the advantage that they are much more symmetric than WH -histograms 
(see Fig. |6|). The histograms were reweighted to the 'equal height' temperature. The 
interface tension can then be obtained from ln[P max /P min ]/(2L^,) — > cr/T, when the 
volume — > oo (cf. Eq. ( |4.5| )). 

In practice, the infinite volume value of a is reached in such large volumes that a 
careful finite size analysis is necessary. Following p4, [32], we fit the data with the 
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Figure 14: The V — > oo extrapolations of the interface tension for mu = 50 GeV (top) 
and 65 GeV (bottom). The values shown have been transformed with Eq. (|7.2j) . 



ansatz 
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^{L x g"g^) 2 -P m in (Lx9s3) 2 



3 11 

- ln(L z g 2 S3 ) - - \n(L x g 2 S3 ) + -G + const. 



(7.2) 
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Figure 15: The infinite volume values for the interface tension. 



The function G interpolates between lattice geometries; the limiting values are G = In 3 
for cubical volumes (L z = L x ) and G = for long cylinders (L z ^> L x ). 

The results of the fit to Eq. ( |7.2| ) are shown in Fig. [TJ] for rhjj = 50 GeV and 65 GeV. 
At 50 GeV the interface tension is weak enough that we are forced to use relatively 
large lattices (L x > 7 g$l) before the ansatz ( [7.2| ) can be used. In contrast, at 65 GeV 
we have an excellent fit to all cylindrical lattices. The small cubical lattices at f3$ = 12 
do not display the flat region in the histograms and are excluded here. 

The V — > oo results are shown in Fig. [15| and in Table [3| The symmetric <-> broken 
U -transition is so weak that our volumes in Table ^| are much too small for a reliable 
interface tension estimate, and we did not attempt it here. In addition, at rhjj = 70 GeV 
at the broken U <-» broken H transition our lattices are too small for a reliable estimate 
of a — in this case the problem is simply that the transition is extremely strong, and 
the tunnelling times are very long even with the multicanonical algorithm. The largest 
lattice we used at my = 70 GeV is 16 3 , and the result for o is to be understood only 
as a qualitative estimate. 
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rhu/GeV 




50 


60 


65 


67 


T c /GeV 


Ps = 12 


90.0073(65) 


87.347(10) 


85.614(19) 


84.538(23) 




/3s = 20 


89.787(16) 




85.324(44) 






continuum 


89.457(41) 




84.89(11) 






perturb. 


92.68 


90.20 


88.56 


87.75* 


vh/T c 


ft? = 12 


1.5757(49) 


1.7838(47) 


1.9113(57) 


2.011(15) 




/3s = 20 


1.5343(89) 




1.892(16) 






continuum 


1.480(23) 




1.864(41) 






perturb.** 


1.327 


1.517 


1.647 


1.715* 


L/T* 


/5s = 12 


1.2477(62) 


1.5499(88) 


1.7357(34) 


1.692(58) 




/3s = 20 


1.2073(91) 




1.692(20) 






continuum 


1.147(25) 




1.626(50) 






perturb. 


0.794 


1.024 


1.178 


1.238* 




/3s = 12 


0.1093(25) 


0.1662(51) 


0.2447(42) 


0.301(14) 




/3s = 20 


0.1043(74) 




0.2346(96) 






continuum 


0.097(19) 




0.219(25) 






perturb. 


0.061 


0.105 


0.151 


0.183* 






symm. <-> broken [7 


broken U <- 


y broken H 


rhu/GeV 




68 


70 


68 


70 


T c /GeV 


/3s = 12 


85.730(23) 


88.471(49) 


81.773(62) 


73.65(19) 




perturb. 


94.90 


98.07 


76.90 


70.05 


v H /T c 


/3s = 12 






2.294(39) 


2.914(76) 




perturb.** 






2.67 


3.16 


vu/T c 


/3s = 12 


1.3415(68) 


1.3663(72) 


1.816(16) 


2.372(31) 




perturb.** 


1.31 


1.30 


2.16 


2.57 




/3s = 12 


0.560(17) 


0.571(16) 


0.9570(58) 


1.402(49) 




perturb. 


0.668 


0.663 


1.434 


2.041 


a/T? 


Ps = 12 






0.877(37) 


1.426(51) 




perturb.*** 






2.9 


5.3 



Table 3: The infinite volume and continuum extrapolations. The continuum values 
have been linearly extrapolated from /3s = 12, 20. The perturbative values for rhjj = 67 
GeV (*) correspond to a transition deep in the (perturbative) metastability region, see 
Fig. |9|. As explained in the text, it should be noted that the perturbative definitions 
for vh /T c ,vu /T c (**) are not exactly the same as the non-perturbative ones. The 
perturbative values (***) for cr/T c 3 represent an upper bound as explained in [T7 . 
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7.5 Correlation lengths 



We measure the screening masses of the operators in Eqs. ( |4.8|) from an additional 
series of runs at rhy = 60 and 68 GeV, using (3$ = 12, 32 2 x 64 lattices. The correlation 
functions are measured in the direction of the X3-axis, and in order to enhance the 
projection to the ground states, we use recursive gauge invariant blocking of the fields 



in the (xi, X2)-plane. The blocking we use is similar to the one in Ref. |[38|| . The fields 
on the level (k + 1) are effectively defined only on the even points of the (fc)-level 
lattice on the (x, y)-plane. The blocking is performed with the transformations (here 
<f) is either H or U, and U = U w or U s , correspondingly) 

<P (k+1 \y) = U {k \x) + \ £ u\ k \x)^\x + i\ 

5 i=±l,2 

U? +1 \y) = U^(x)U^ k \x + i), (7.3) 

U: {k \x) = l -U\ k \x)+ l - J2 uf\x)U?\x + j)uf ) \x + i) t 

j=±i' 

where (xi, x%, x$) = (2yi, 2y 2 , 2/3) and i — 1, 2, z' = 3 — i. We use the blocked fields to 
construct the blocked SU(2) operators 



x 



scalar: S { h\x) = H^{x)H {k) { 

vector: V^ k \x) = luiH^ (x)UY {k \x)H^ k \x + i), (7.4) 

++ glueball: G%\x) = 1 - ±Tr P^ (fc) (x), 

where is the (xi, X2)-plane plaquette (resp. for SU(3)). The operators are summed 
over (xi, X2)-planes at each value of x%, and we measure the plane-plane correlation 
functions. The operators are blocked up to 5 times, and we measure the correlation 
functions for each level separately. The masses are read from the exponential fall-off of 
these functions. Due to the periodicity in the X3 direction, we fit a hyperbolic cosine to 
the vector channel and a constant + hyperbolic cosine to the scalar channel correlation 
functions. All of the fits use the full covariance matrix of the correlation functions. 
The fitting range is automatically selected so that the range is as long as possible while 
still keeping the confidence level acceptable. 



The screening masses for rrtu = 60 GeV are shown in Fig. |16| and for mu = 68 GeV in 



Fig. [17|. For each temperature and operator we choose the blocking level which yields 
the best result. The optimal level was usually 4 or 5 (level 1 meaning no blocking). 
We observe the following: 

(1) The scalar operators Sh, Sjj, Gh and G\j all have the same quantum numbers 
and couple to the same set of states. A reliable measurement of the scalar mass 
spectrum would require the diagonalization of the full cross-correlation matrix of the 



scalar states, as was done for SU(2)+Higgs by Philipsen et al. [[35J. Lacking the cross- 



correlation matrix, we were not able to reliably extract the glueball masses. Moreover, 
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Figure 16: The screening masses at ihjj = 60GeV, measured from a (3s = 12, 32 2 
lattice. The transition temperature is marked with a vertical dashed line. 



x 64 



at fhu = 60 GeV both the Sh and Su operators give the mass of the lightest scalar 
state Si/-Q 

(2) The vector operators do not couple to each other in a similar way, and we can 
measure the corresponding masses in all cases. 

(3) The masses show a clear discontinuity at the critical temperatures. For fhu = 
60 GeV the behaviour of Sh and V H at the transition temperature is strongly rem- 
iniscent of that for the corresponding SU(2)+Higgs theory masses |32], [35| . At the 
same time, the mass of the SU(3) vector state V v increases dramatically when the 
system undergoes the transition from the symmetric phase into the broken H phase. 
This can be qualitatively understood in terms of confinement and bound states (see 
also |55|, |56"1|): when the if-field enters the broken phase, the effective mass term of the 



[/-field increases substantially, due to the jH^HU^U -coupling. Thus the mass of the 
'heavy squark' meson-like bound state Wlf also increases. 



8 The state is naturally a mixture of all scalar operators, but since the dominant coupling is to the 
Sh = H^H operator, we keep this name for the state. 
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Figure 17: The screening masses at my = 68GeV, measured from a (3$ — 12, 32 2 x 64 
lattice. The transition temperatures are marked with the vertical dashed lines. 
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8 Discussion of the results 



Let us then discuss the comparison of the lattice results with perturbation theory. 

The phase diagram and the critical temperatures are shown in Fig. ||]. It is seen 
that the phase diagram is qualitatively the same as in perturbation theory, although 
the critical temperatures and the triple point have been displaced by a few GeV. The 
displacement is statistically significant, as the errors of the lattice points are very small 
even after continuum extrapolation. We have no clear theoretical explanation for the 
discrepancy between lattice results and perturbation theory: the reason might be, e.g., 
a three-loop perturbative effect, or a genuine non-perturbative contribution (see, e.g., 
|57| , |58||). Let us note that the non-perturbative critical temperature for the electroweak 
phase transition is smaller than the perturbative one also for the Standard Model . 

The main result of this paper is Fig. [H| which shows the latent heat. The latent 
heat is the most important gauge invariant physical characterization of the strength 
of a first order transition. We observe that the non-perturbative transition to the 
standard electroweak minimum at my < 67 GeV is significantly (up to 45%) stronger 
than the perturbative transition. This behaviour is in strong contrast to that in the 
Standard Model, where the latent heat is smaller than in perturbation theory at least 
for Higgs masses m# > 70 GeV J32, |33|. Again, we have no clear explanation for this 
behaviour. 

In the regime of the two-stage transition (my > 67 GeV), the comparison with per- 
turbation theory is not quite as straightforward, as the whole pattern has been shifted 
in my. Nevertheless, the first stage of the transition (symmetric — > broken U) is clearly 
weaker than in perturbation theory, analogously to what happens for the normal tran- 
sition in the Standard Model. 

A similar pattern as for the latent heat, is observed for the interface tension, Fig. |L5 . 
The relative strengthening effect at my < 67 GeV is of the same order of magnitude 
as for the latent heat. In the regime of large my, an important observation is that 
even non-perturbatively the interface tension grows very fast with increasing my. The 
qualitative similarity with the perturbative estimate is impressive, taking into account 



that the perturbative numbers shown represent a very rough upper bound estimate [17 



The scalar field expectation values are shown in Fig. O. To get a transition strong 
enough for baryogenesis, one needs to have Vh/T >1[|. Again, we observe a value 
larger than in perturbation theory in the regime rhy < 67 GeV, and a rapid increase in 
vh/T c in qualitative accordance with perturbation theory in the regime of the two-stage 
transition, my > 67 GeV. It should be kept in mind, though, that the definitions of the 
objects shown are not exactly the same in perturbation theory and on the lattice, as 
has been discussed in Sec. £|. 

Finally, consider the correlation lengths in Figs. [HJ |T7]. For the correlation lengths, 

9 To be more precise, the requirement should probably be vh/T> 1.2 . . . 1.5 |52|. A recent lattice 
computation B9[ favours the lower end of this range. 
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the non-perturbative confining dynamics of the symmetric phase shows up in a very 
dramatic way, and comparison with perturbation theory is only possible in the broken 
phases. Thus, for T < T c 2 in Fig. [17|, one can compare the Higgs and SU(2) vector 
masses with perturbation theory, and for T c 2 < T < T c l , the stop and SU(3) vector 
masses. In these regimes, we observe that compared with the tree-level perturbative 
masses (m H /g 2 S3 ~ 1.2, m w /g 2 S3 ~ 0.8 for T < T c , 2 , andm iR /g 2 S3 ~ 1.4 . . . 0.9, m G /g 2 S3 ~ 
1.2 .. . 0.7 for T C} 2 < T < T c i), the non-perturbative scalar masses are somewhat smaller 
and the non-perturbative vector masses somewhat larger. This is in accordance with 
the pattern observed for the Standard Model in p2| |. In the other regimes where the 



excitations feel an unbroken gauge group, the physical masses are those of bound states 
and very large. 



9 Conclusions 

We have studied the electroweak phase transition in the MSSM with non-perturbative 
lattice Monte Carlo simulations, in the regime of large (m# « 95 GeV) Higgs masses 
and small {mi R ~ 150 . . . 160 GeV) stop masses. Several properties of the phase tran- 
sition have been determined: the phase structure and critical temperatures, the latent 
heat, the interface tension, and the correlation lengths. The results have been extra- 
polated to the infinite volume and continuum limits. 

The main conclusion of the study is that at least for the parameter values used, the 
electroweak phase transition in the MSSM is significantly stronger than indicated by 
2-loop Landau-gauge ft = T perturbation theory. If the same pattern remains there 
for larger Higgs masses, then this implies that previous perturbative Higgs and stop 
mass bounds for electroweak baryogenesis are conservative estimates. In particular, 
the electroweak phase transition would then be strong enough for baryogenesis for all 
allowed Higgs masses in this regime (m# < 105 GeV) UlSR . 

This result certainly provides a strong additional motivation for experimental Higgs 
and stop searches at LEP and the Tevatron ]Tj|. Moreover, it provides a strong moti- 



vation for more precise studies of the non-equilibrium CP- violating real time dynamics 
and baryon number generation at the electroweak phase transition. It would also be 
interesting (and straightforward) to extend the present simulations to other parame- 
ter values, such as a Higgs mass very close to the upper bound mu ~ 105 GeV, and 
non-vanishing mixing in the squark sector. In addition, it would be useful to have 
an explanation for why the non-perturbative transition is stronger than indicated by 
2-loop perturbation theory in the present case, in contrast to the SU(2)+Higgs model. 

For the largest values of rhu-, we have observed that the electroweak phase transition 
can take place in two stages, and we have analyzed this regime in detail. The values 
needed for rhu are somewhat larger than in perturbation theory, corresponding to 
smaller values of mz . The main physical implication of the two-stage transition is that 
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it is a way of making the transition where the Higgs field gets a vacuum expectation 
value and the sphaleron rate is switched off, extremely strong. The intermediate regime 
where the stop field has an expectation value, might also have exotic properties. 

At the same time, the price to be paid for a strong transition is that the interface 
tension is quite large. This implies that the supercooling taking place is significant, of 
order 35% already at rhu = 68 GeV (the nucleation temperature T n can be estimated 
from 1 — T n /T c oc a 3 / 2 / (LT^ 2 ); see, e.g, |T7| ). As the supercooling is getting larger, 



there is the danger that the transition does not take place at all during cosmological 
time scales, which would forbid this scenario. Thus the results can also be interpreted 
as an upper bound for rhu, or a lower bound for m f - fl 



17 . Note, however, that non-zero 



squark mixing parameters seem to significantly reduce the possibility of the two-stage 
transition and at the same time they allow for smaller stop masses. 



Finally, let us note that the latent heat and interface tension determined here can 
also be used for estimates of the non-equilibrium real-time dynamics of the transition, 
such as the nucleation temperature (see above), the velocities of expanding bubbles, 
whether reheating to T c takes place after the transition, etc [30|. For instance, using 
the non-perturbative values of L/T^ and cr/T c 3 for rhu — 50 GeV, it appears that 
reheating to T c does not take place, so that the scalar vacuum expectation value after 
the transition is even larger than that at T c 
suppress the sphaleron rate. 



given in Fig. O. This serves to further 



Acknowledgements 

The simulations were made with a Cray T3E at the Center for Scientific Comput- 
ing, Finland. We acknowledge useful discussions with K. Kainulainen, K. Kajantie, 
M. Losada, CD. Moore, T. Prokopec, M. Shaposhnikov and C. Wagner. This work 
was partly supported by the TMR network Finite Temperature Phase Transitions in 
Particle Physics, EU contract no. FMRX-CT97-0122. 



References 

[1] V.A. Kuzmin, V.A. Rubakov and M.E. Shaposhnikov, Phys. Lett. B 155 (1985) 
36; M.E. Shaposhnikov, Nucl. Phys. B 287 (1987) 757. 



[2] V.A. Rubakov and M.E. Shaposhnikov, Usp. Fiz. Nauk 166 (1996) 493 jhcp- 



ph/9603208 



[3] K. Kajantie, M. Laine, K. Rummukainen and M. Shaposhnikov, Phys. Rev. Lett. 
77 (1996) 2887. 

[4] M. Gurtler, E.-M. Ilgenfritz and A. Schiller, Phys. Rev. D 56 (1997) 3888. 



38 



[5] M. Giovannini and M.E. Shaposhnikov, Phys. Rev. D 57 (1998) 2186. 

[6] S. Myint, Phys. Lett. B 287 (1992) 325; G.F. Giudice, Phys. Rev. D 45 (1992) 
3177. 

[7] J.R. Espinosa, M. Quiros and F. Zwirner, Phys. Lett. B 307 (1993) 106. 

A. Brignole, J.R. Espinosa, M. Quiros and F. Zwirner, Phys. Lett. B 324 (1994) 
181. 

[9] M. Carena, M. Quiros and C.E.M. Wagner, Phys. Lett. B 380 (1996) 81. 
[10] J.R. Espinosa, Nucl. Phys. B 475 (1996) 273. 

[11] D. Delepine, J.-M. Gerard, R. Gonzalez Felipe and J. Weyers, Phys. Lett. B 386 
(1996) 183. 



[12] M. Laine, Nucl. Phys. B 481 (1996) 43 ||hep-ph/9605"283| 



[13] J.M. Cline and K. Kainulainen, Nucl. Phys. B 482 (1996) 73; Nucl. Phys. B 510 
(1998) 88. 

[14] M. Losada, Phys. Rev. D 56 (1997) 2893; G.R. Farrar and M. Losada, Phys. Lett. 
B 406 (1997) 60. 

[15] J.M. Moreno, D.H. Oaknin and M. Quiros, Nucl. Phys. B 483 (1997) 267; Phys. 
Lett. B 395 (1997) 234. 

[16] B. de Carlos and J.R. Espinosa, Nucl. Phys. B 503 (1997) 24. 

[17] D. B6deker, P. John, M. Laine and M.G. Schmidt, Nucl. Phys. B 497 (1997) 387 
||hep-ph/96l2364H . 



[18] M. Carena, M. Quiros and C.E.M. Wagner, CERN-TH/97-190 ||hep-ph/ 9 7102011 



[19] M. Carena, M. Quiros, A. Riotto, I. Vilja and C.E.M. Wagner, Nucl. Phys. B 503 
(1997) 387; T. Multamaki and I. Vilja, Phys. Lett. B 411 (1997) 301. 

[20] J.M. Cline, M. Joyce and K. Kainulainen, Phys. Lett. B 417 (1998) 79. 



[21] K. Enqvist, A. Riotto and I. Vilja, OUTP-97-49-P ||hep-ph/9710373l ; A. Riotto 



CERN-TH/98-81 ||hep-ph/9803357ll , and references therein. 



[22] M.P Worah, Phys. Rev. D 56 (1997) 2010; Phys. Rev. Lett. 79 (1997) 3810. 
[23] J.M. Moreno, M. Quiros and M. Seco, IEM-FT-168-98 |[hep-ph/ 980 1272] . 



39 



[24] K. Funakubo, A. Kakuto, S. Otsuki and F. Toyoda, SAGA-HE-132 jhep- 
| ph/9802276[ 



M. Laine and K. Rummukainen, CERN-TH/98-121 ||hep-ph/9804255[ . 

F. Csikor, Z. Fodor, J. Hein, K. Jansen, A. Jaster and I. Montvay, Phys. Lett. 
B 334 (1994) 405; Z. Fodor, J. Hein, K. Jansen, A. Jaster and I. Montvay, Nucl. 
Phys. B 439 (1995) 147; F. Csikor, Z. Fodor, J. Hein and J. Heitger, Phys. Lett. 
B 357 (1995) 156; F. Csikor, Z. Fodor, J. Hein, A. Jaster and I. Montvay, Nucl. 
Phys. B 474 (1996) 421; Y. Aoki, Phys. Rev. D 56 (1997) 3860. 

P. Ginsparg, Nucl. Phys. B 170 (1980) 388; T. Appelquist and R. Pisarski, Phys. 
Rev. D 23 (1981) 2305; S. Nadkarni, Phys. Rev. D 27 (1983) 917. 

K. Farakos, K. Kajantie, K. Rummukainen and M. Shaposhnikov, Nucl. Phys. 
B 425 (1994) 67; K. Kajantie, M. Laine, K. Rummukainen and M. Shaposh- 
nikov, Nucl. Phys. B 458 (1996) 90 |hep-ph/95Qg57S|| ; Phys. Lett. B, in press 
|hep-ph/971053qi . 



A. Jakovac, K. Kajantie and A. Patkos, Phys. Rev. D 49 (1994) 6810; A. Jakovac 
and A. Patkos, Phys. Lett. B 334 (1994) 391; Nucl. Phys. B 494 (1997) 54. 

E. Braaten and A. Nieto, Phys. Rev. D 51 (1995) 6990; Phys. Rev. D 53 (1996) 
3421. 

K. Kajantie, K. Rummukainen and M. Shaposhnikov, Nucl. Phys. B 407 (1993) 
356. 

K. Kajantie, M. Laine, K. Rummukainen and M. Shaposhnikov, Nucl. Phys. B 
466 (1996) 189. 

E. -M. Ilgenfritz, J. Kripfganz, H. Perlt and A. Schiller, Phys. Lett. B 356 (1995) 
561; M. Giirtler, E.-M. Ilgenfritz, J. Kripfganz, H. Perlt and A. Schiller, Nucl. 
Phys. B (Proc. Suppl.) 49 (1996) 312 ||hep-lat/9512022| ; Nucl. Phys. B 483 (1997) 
383; M. Giirtler, E.-M. Ilgenfritz and A. Schiller, Eur. Phys. J. C 1 (1998) 363. 

F. Karsch, T. Neuhaus, A. Patkos and J. Rank, Nucl. Phys. B 474 (1996) 217. 

O. Philipsen, M. Teper and H. Wittig, Nucl. Phys. B 469 (1996) 445; OUTP-97- 
44-P |iep-lat/9709145| . 

CD. Moore and N. Turok, Phys. Rev. D 55 (1997) 6538. 
M. Laine and O. Philipsen, Nucl. Phys. B, in press ||hep-lat/9711022 



40 



[38] K. Kajantie, M. Laine, K. Rummukainen and M. Shaposhnikov, Nucl. Phys. B 
493 (1997) 413 [ |hep-lat/ 961 2006; ]. 

[39] T. Reisz, Z. Phys. C 53 (1992) 169; P. Lacock, D.E. Miller and T. Reisz, Nucl. 
Phys. B 369 (1992) 501; L. Karkkainen, P. Lacock, D.E. Miller, B. Petersson and 
T. Reisz, Phys. Lett. B 282 (1992) 121; Nucl. Phys. B 418 (1994) 3; L. Karkkainen, 
P. Lacock, B. Petersson and T. Reisz, Nucl. Phys. B 395 (1993) 733. 

[40] A. Hart, O. Philipsen, J.D. Stack and M. Teper, Phys. Lett. B 396 (1997) 217. 

[41] K. Kajantie, M. Laine, K. Rummukainen and M. Shaposhnikov, Nucl. Phys. B 
503 (1997) 357 ||hep-ph/97044T6| ; K. Kajantie, M. Laine, J. Peisa, A. Rajantie, 
K. Rummukainen and M. Shaposhnikov, Phys. Rev. Lett. 79 (1997) 3130 ||hep-| 

| ph/9708207| . 

[42] M. Karjalainen and J. Peisa, Z. Phys. C 76 (1997) 319 ||hep-lat/9607023l ; K. 
Kajantie, M. Karjalainen, M. Laine and J. Peisa, Phys. Rev. B 57 (1998) 3011 
|cond-mat/970405q] ; Nucl. Phys. B, in press ||hep-lat/97lT048|] ; K. Kajantie, M. 



Karjalainen, M. Laine, J. Peisa and A. Rajantie, Phys. Lett. B, in press ||hep- 
ph/9803367|| . 



[43] P. Dimopoulos, K. Farakos and G. Koutsoumbas, Eur. Phys. J. C 1 (1998) 711. 
[44] M. Losada, to appear. 

[45] S. Dimopoulos, S. Raby and L. Susskind, Nucl. Phys. B 173 (1980) 208. 

[46] K. Farakos, K. Kajantie, K. Rummukainen, and M. Shaposhnikov, Nucl. Phys. B 
442 (1995) 317 ||hep-lat/941209ll ; M. Laine, Nucl. Phys. B 451 (1995) 484 jlep 
| lat/9504Q0l|| . 



[47] M. Laine and A. Rajantie, Nucl. Phys. B 513 (1998) 471 |hep-lat/9705D03fl . 
[48] G.D. Moore, McGill-97-23 [frep-lat/97090531 . 

[49] G.D. Moore, Nucl. Phys. B 493 (1997) 439; Phys. Lett. B 412 (1997) 359. 
[50] A.D. Kennedy and B.J. Pendleton, Phys. Lett. B 156 (1985) 393. 
[51] B. Berg, J. Stat. Phys. 82 (1996) 331. 

[52] The MILC collaboration public domain QCD code is available at 
Qhttp : //physics . indiana. edu/~sg/milc .htmlp . 

[53] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635. 



41 



[54] Y. Iwasaki, K. Kanaya, L. Karkkainen, K. Rummukainen and T. Yoshie, Phys. 
Rev. D 49 (1994) 3540. 

[55] H.-G. Dosch, J. Kripfganz, A. Laser and M.G. Schmidt, Phys. Lett. B 365 (1996) 
213; Nucl. Phys. B 507 (1997) 519. 

[56] W. Buchmiiller and O. Philipsen, Phys. Lett. B 397 (1997) 112. 

[57] M. Shaposhnikov, Phys. Lett. B 316 (1993) 112. 

[58] S.J. Huber, A. Laser, M. Reuter and M.G. Schmidt, to appear. 

[59] G.D. Moore, McGill-97-36 |frep-ph/9801204 |. 

[60] H. Kurki-Suonio and M. Laine, Phys. Rev. Lett. 77 (1996) 3951. 



42 



